#include "pseudo.hpp"
#include <iostream>

typedef int (* SELECT)(std::complex<double> *);

extern "C" void 
zgees_(char *jobvs, char *sort, SELECT select, int *n,
       DComplex *a, int *lda, int *sdim, DComplex *w, DComplex *vs, int *ldvs,
       DComplex *work, int *lwork, double *rwork, int *bwork, int *info);

void 
zgees(char jobvs, char sort, SELECT select, int n, DComplex *a, int lda,
      int &sdim, DComplex *w, DComplex *vs, int ldvs, DComplex *work,
      int lwork, double *rwork, int *bwork, int &info)
{
  
  zgees_(&jobvs, &sort, select, &n, a, &lda, &sdim, w, vs, &ldvs, work, &lwork,
	rwork, bwork, &info);
}


Pseudospectra::Pseudospectra(DoubleMatrix a):
  _N(a.dim1()), _T(_N, _N), _eig(_N)
{

  init(a);
}

Pseudospectra::Pseudospectra(double *a, int n):
  _N(n), _T(_N, _N), _eig(_N)
{
  DoubleMatrix am(_N, _N, n);
  init(am);
}

void
Pseudospectra::init(DoubleMatrix &a)
{
  for (int i = 0; i < _N; ++i)
    for (int k = 0; k < _N; ++k)
      _T[i][k] = a[k][i];  	// Lapack uses fortran indices.
  
  // Perform Schur decomposition of the matrix:
  int info;
  int sdim;
  int lwork = 4*_N;
  ComplexArray work(lwork);
  DoubleArray rwork(_N);
  
  zgees('N', 'N', 0, _N, _T[0], _N, sdim, _eig, 0, 1, work, lwork, rwork, 0, info);


  std::cout << _T;
}
  
double
Pseudospectra::sigmin(double x, double y, ComplexMatrix T1)
{

  // The correct matrix is 
}

int main(){

  /*
  A =

   0.710806  -0.031361   0.412391  -0.120999
  -1.580280  -0.838259   0.807765  -0.649884
   0.256918   0.613558  -0.029110  -0.859698
   0.521888  -0.794682  -0.462181  -0.325377
  */

  DoubleMatrix A(4,4);
  A[0][0] =  0.710806;
  A[0][1] = -0.031361;
  A[0][2] = 0.412391;
  A[0][3] = -0.120999;
  A[1][0] =  -1.580280;
  A[1][1] = -0.838259;
  A[1][2] = 0.807765;
  A[1][3] = -0.649884;
  A[2][0] =  0.256918;
  A[2][1] = 0.613558;
  A[2][2] = -0.029110;
  A[2][3] = -0.859698;
  A[3][0] = 0.521888;
  A[3][1] = -0.794682;
  A[3][2] = -0.462181;
  A[3][3] = -0.325377;
  


  Pseudospectra pseudo(A);
  
}


